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Abstract 

Control-based continuation is technique for tracking the solutions and bifurcations of nonlinear experiments. 
The basic idea is to apply the method of numerical continuation to a feedback-controlled physical experiment. 
Since in an experiment it is not (generally) possible to set the state of the system directly, the control target is 
used as a proxy for the state. The challenge then becomes to determine a control target such that the control 
is non-invasive, that is, it stabilises a steady-state (or periodic orbit) of the original open-loop experiment 
without altering it otherwise. 

Once implemented, control-based continuation enables the systematic investigation of the bifurcation 
structure of a physical system, much like if it was numerical model. However, stability information (and hence 
bifurcation detection and classification) is not readily available due to the presence of stabilising feedback 
control. This paper uses methods from the system identification community to extract stability information 
in the form of Floquet multipliers from the closed-loop experiment, thus enabling the direct detection of 
bifurcations. In particular, it is shown that a periodic auto-regressive model with exogenous inputs (ARX) 
can be constructed that approximates the time-varying linearisation of the experiment around a particular 
periodic orbit. This method is demonstrated using a physical nonlinear tuned mass damper. 

Keywords: numerical continuation, bifurcation theory, system identification, feedback control 


1. Introduction 

Control-based continuation is a systematic method for performing bifurcation studies on physical ex¬ 
periments. Based on modern feedback control schemes it enables dynamical phenomena to be detected 
and tracked as system parameters are varied in a similar manner to how nonlinear numerical models can 
be investigated using numerical continuation. Control-based continuation was originally developed as an 
extension to Pyragas’ time-delayed feedback control [1-3] to make it more robust and suitable for parameter 
studies [4], though its current implementation contains no elements of time-delayed feedback. 

The use of feedback control for the investigation of nonlinear systems is not new; in addition to time- 
delayed feedback, methods such as OGY control [5] have previously been employed to provide non-invasive 
control to stabilise unstable orbits and investigate dynamical phenomena. Indeed, previous authors have 
gone as far as to implement such control schemes within parameter continuation studies to numerically 
simulate particular experiments such as atomic force microscopes [6] or reaction kinetics [7]. Control-based 
continuation goes beyond these particular methods to allow the use of almost any feedback control scheme 
and, as such, it is a general purpose tool applicable to a wide range of physical experiments. 

Control-based continuation has been successfully applied to a range of experiments including a paramet¬ 
rically excited pendulum [8], nonlinear energy harvesters [9, 10] and a bilinear oscillator ]11, 12]. In each 
case, periodic orbits have been tracked through instabilities such as saddle-node bifurcations (folds) thus 
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revealing a great deal of dynamical information about the system in question, including the location of a 
codimension-2 cusp bifurcation in one case [13]. As well as bifurcations, other dynamic features such as 
backbone curves can also be tracked with control-based continuation [14]. 

Although the basic scheme for control-based continuation is well established (an overview is provided in 
Sec. 2), it lacks many of the features of standard numerical continuation schemes such as bifurcation detection. 
Only saddle-node bifurcations (folds) can be detected readily and that is because they are geometric features 
in the solution surface. Bifurcations such as period-doubling bifurcations are not geometric features and can 
go undetected due to the stabilising affect of the feedback controller. Similarly, the inclusion of the feedback 
controller means that methods for calculating eigenvalues/Floquet multipliers and basins of attraction from 
experiments such as [15-17] are not helpful since they indicate the stability of the closed-loop system rather 
than the open-loop system. 

In this paper we consider only periodically forced systems and hence study periodic orbits, though 
there is no reason that the methods developed should not be applicable to autonomous systems as well. 
In Sec. 3, we propose a method for calculating the stability (the Floquet multipliers and associated stable 
and unstable eigendirections) based on the estimation of a local linearisation around a stabilised periodic 
orbit. We demonstrate the effectiveness of this method in Sec. 4 and Sec. 5 by applying it to a (physical) 
nonlinear mass-spring-damper system where the nonlinearity is geometric in nature — the springs are 
mounted perpendicular to the direction of motion. 

2. Control-based continuation 

Numerical continuation is a path following method used to track solution branches as parameters of the 
system in question are varied. In a nonlinear system, these solution branches can encounter bifurcations at 
particular parameter values which results in a qualitative change in the dynamics of the system. Numerical 
continuation enables these bifurcations to be detected and tracked in turn. It is typically applied to differential 
equation models but it can be used more widely, for example on finite element models. 

At a basic level, numerical continuation tracks the solutions of an arbitrary nonlinear function, a zero 
problem given by 

fix, A) = 0, / : K” X MP ^ (1) 

where x is the system state and A is the system parameter(s). A common example of this is tracking the 
equilibria of an ordinary differential equation with respect to a single parameter. In this case n = m and 
p = I in Eq. (1), that is, / = 0 defines a one-dimensional curve. Alternatively the function / can arise 
from the discretisation of a periodic orbit. Numerical continuation works in a predictor-corrector fashion; 
at each step a new solution x is predicted from previously determined solutions and then the solution is 
corrected using a nonlinear root finder applied to the function / (typically a Newton iteration). The use 
of a nonlinear root finder means that the stability or instability of solutions is unimportant. In certain 
circumstances (for example, near a fold or saddle-node bifurcation) the function / must be augmented with an 
additional equation — the pseudo-arclength equation — which enables the numerical continuation scheme to 
track solution curves that double back on themselves. In these circumstances, without the pseudo-arclength 
equation the correction step for a fixed set of parameter values A will fail since no solution exists. For 
extensive information and guidance on numerical continuation see the textbooks [18, 19]. Numerical software 
is also readily available in the form of CoCo [20] and AUTO [21] amongst others. 

Control-based continuation is a means for defining a zero-problem based on the outputs of a physical 
experiment, thus enabling numerical continuation to be applied directly without the need for a mathematical 
model. To do this there are two key challenges to overcome. 1) In general, it is not possible to set the state 
X of the physical system and so it is not possible to evaluate / at arbitrary points. 2) The physical system 
must remain around a stable operating point while the experiment is running. While a numerical model 
going unstable might prove to be a mild annoyance, a physical system going unstable can prove dangerous. 

In order to overcome these challenges, a feedback controller is used to stabilise the system and the control 
target (or reference signal) acts as a proxy for the system state. The feedback controller takes the form 


uit) = gix*{t) - x{t)) 
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( 2 ) 


where x*{t) is the control target and g is a suitable control law such as proportional-derivative (PD) control 
(as used in this paper) where 

u{t) = Kp{x*{t) - x{t)) -I- Kd{x*{t) - x{t)). (3) 

For the method outlined in this paper, the choice of control law is at the discretion of the user; any suitable 
stabilising feedback control scheme can be used. The challenge here is to devise a scheme for embedding the 
feedback control within the numerical continuation such that the controller becomes non-invasive, that is, the 
controller does not affect the locations of any invariant sets in the experiment such as equilibria or period 
orbits. This requirement for non-invasiveness defines the zero problem; a control target must be chosen such 
that the control action 

u{t) = 0. (4) 

In this paper, we consider the case of a periodically forced experiment with forcing frequency w and, as 
such, only consider periodic motions. In this case it is appropriate to consider a Fourier discretisation of Eq. (4) 
and so find the coefficients of the Fourier series of the control target a:*(t) = /2 + Aj cos{ju!t) + 

Bj sm{jujt) such that Eq. (4) is satisfied. (In other circumstances different discretisations may be appropriate.) 
In this case the control action u has a Fourier series representation given by 

AU ™ 

u{t) = -^ + X! cos(ja;f) -f BJ sin{jujt). (5) 

i=i 

where the Fourier coefficients AJ and BJ are derived directly from the measured control action Eq. (2), that 
is, 

27r 

AJ = — f g{x*{t) — x{t)) cos{jLot) dt, 
tt Jo 

2-7r 

BJ = - [ g{x*{t) - x(t)) sin(jwt) dt, 

Jo 

Hence the discretised zero problem is defined as 

0 = AJ, 0 = BJ Vj. 


for j = 0,1,2,... 

(6a) 

for j = 1,2,.... 

(6b) 


( 7 ) 


To solve Eq. (7) standard root finding algorithms can be used, however, any required derivatives must be 
estimated using finite differences from experimental data after adjusting the experiment inputs. Consequently, 
gradient-based methods can be slow despite their good convergence rates. In previous publications a 
Newton-Broyden method, which avoids recomputing derivative information for successive iterates, has proven 
effective [8, 10, 11]. 

In this paper, since the control acts through the same mechanism as the forcing, we are able to use a 
quicker method which exploits the fact that we are performing a parameter study in the forcing amplitude [13]. 
Consider the case where the total input to the system is given by 

i(t) = pit) + u{t) (8) 

where pit) is the forcing signal and u(t) is the control action. Furthermore, we consider the case of sinusoidal 
forcing where 

pit) = acos(wt) -I- &sin(tLit). (9) 

For an arbitrary control target u*it) and forcing input pit), the Fourier coefficients of Eq. (5) will be 
non-zero. However, the contribution in the fundamental mode (coefficients AJ and B'^) can be lumped 
together with the coefficients a and b of the forcing term giving a new effective forcing amplitude of 

F = ^Jia + Aj)"^ + ib + BJ)"^. 
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Hence, once the higher Fourier modes of the control action u{t) are set to zero (as described below), the total 
input to the system will be i{t) = T cos{ujt + (f>) (the phase (j) is unimportant since the system is time invariant). 
In essence, instead of setting the forcing amplitude and trying to calculate the correct corresponding control 
target, we set the control target and measure the corresponding forcing amplitude. 

Though this procedure leaves the Aq and the higher harmonics untouched, the corresponding control 
target coefficients required to set the control action to zero can be quickly determined using a fixed-point 
iteration. In an iterative manner the remaining coefficients of the control target x*(t) are simply set equal to 
the measured coefficients of the system response x{t), the iteration finishes when the response and the control 
target remain equal for a certain period of time. For the system described below, this takes a single iteration. 

A fixed-point iteration cannot be applied to the coefficients of the fundamental mode since, generically, 
instabilities in the system will manifest in the fundamental mode. 

It is important to emphasise that this procedure does not depend on the specifics of the control law used 
in the experiment. All that is required is that smooth changes in the control target x*{t) result in smooth 
changes in the Fourier coefficients Eq. (6). As such, this method is convenient in realistic settings where the 
control law is more sophisticated than simple PD control and filtering of the signal is required; these effects 
are simply captured in the control action u{t) and fed into Eq. (6). 


3. Identification of a linearisation 

For a typical ordinary differential equation model, the right-hand side of which is given by h(x(t)), the 
Floquet multipliers and hence stability of a periodic orbit are determined by integrating over one period the 
first variational equation 

^ = A{x{t))y{t) (11) 

where y{t) represents the deviation from a predetermined periodic orbit x{t) and the Jacobian matrix A is 
calculated from the derivatives of h with respect to x evaluated along the periodic orbit given by x{t) (as 
such A is a time varying quantity). 

In the context of control-based continuation, even determining whether an orbit is stable or unstable is 
problematic due to the presence of stabilising feedback control. In [12] a number of measures are suggested to 
overcome this problem but all require turning off the feedback control for a period of time; in many situations 
this is not desirable as damage could be caused to the experiment or even the experimenter. As such, here 
we attempt to fit a time-varying linearisation using techniques from the system identification community. 
From the fitted time-varying linearisation we are able to calculate the corresponding Floquet multipliers and 
hence the stability properties of the periodic orbit. 

We start by assuming that the experiment of interest is undergoing a periodic motion x(t) for a given 
forcing input i{t) = p{t) + u{t) (cf. Eq. (8)). In order to generate data with which to fit a time-varying 
linearisation, the system is perturbed using filtered Gaussian white noise rj{t) such that the total input to 
the experiment is 

i{t) = i{t) + ^{t) + u{t) - u{t). (12) 

The additional u{t) — u{t) term arises due to the presence of the feedback controller acting against the applied 
perturbation. (Details of rj{t) are below.) Finally, we define the perturbed system response 

y{t) = x{t) - x(t) (13) 

and, similarly, the perturbed system input 

k{t) = i{t) — i{t) . (14) 


Rather than try to fit a continuous time model of the response to perturbations such as Eq. (11), which 
requires the estimation of derivatives from experimental data, we instead fit the coefficients of a discrete-time 
multiple-input multiple-output (MIMO) auto-regressive model with exogenous inputs (ARX) of the form 

B{q-^)y{T) = A(g-i)k(T) + e(T) (15) 
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Figure 1: An example discretisation of the perturbed response of a periodic orbit with m = 5 and n = 2; the period of the orbit 
is normalised to 1. Here y = [j/o^ ^ 2 ? 2 / 3 ? ^^ 4 ]- Since n = 2, yo is purely a function of [ 2 / 1 , 2 / 2 ] and [A:o,^i,fc 2 l- Similarly, 2/4 is 

purely a function of [yoq~^ ,yiq~^] and [k 4 , koq~^^ kiq~^], where q~^ is the backward shift operator. Thus the linearised model 
Eq. (15) allows the construction of a linear period map from yq~^ to y. 


where y{T) = [y{T - i/m)]i^o,,,jn-i, k(T) = [fc(T - and e(T) = [e{T - are 

vectors of data points of the perturbed system response, the perturbed system input (due to the control action 
and an additional random perturbation) and the model error respectively, sampled across a single period of 
oscillation and q~^ is the backward shift operator [22, Sec. 6.2]. Thus y(T) corresponds to a discretisation 
of the perturbed system response using m points over the period; this discretisation is illustrated in Fig. 1. 

and A{q~^) are square m x m matrices of polynomials in here we restrict the polynomials to 
being first order (at most). Thus Eq. (15) acts as a period map with all the dynamics of interest encoded in 
B{q~^). (The matrix A{q~^) though required for system identification are not used to infer stability.) 

The model error e(T) is not measured directly but is minimised by the particular system identification 
method used on Eq. (15). 

ARX models are used extensively in the system identification and time series analysis communities. 
Their simplicity has seen them applied to a wide range of topics. They are particularly appropriate in 
situations where discretely sampled data is available as they avoid the need of estimating derivatives. For 
more information see the textbook by Hamilton [23] or one of the many other books on this topic. 

At each point we assume linear observability to state that the next point in the time evolution of the 
linearisation is determined entirely by the previous n points alone. We assume that n < m (increase the 
value of TO as appropriate) to obtain a banded matrix structure for B{q~^) of 


H(g-i) 


1 h.i 

0 1 


bm.iq ^ bm, 2 q 


and for A{q ^) 


A{q-^) 


01,0 Oip 

0 02,0 

am,iq~^ am,2q' 


bl,n-l 

b2,n-2 

bm,nq 

Ol,n-l 

02,n-2 

^m,nQ 


^l,n 

&2,n-l 

0 

Ol,n 

02,n-l 

0 


0 0 

b2,n 0 

0 0 

0 0 

02,n 0 

0 0 


0 

0 

1 


0 

0 


^m,0 


(16) 


(17) 


Thus the system identification procedure must identify the m{2n + 1) coefficients (aij and bij) within these 
matrices to fully identify the linear model. The optimal values for m and n can be estimated using the Akaike 
information criterion (AIC) [24] or similar. Note that while increasing n increases the data requirements for 
successful system identihcation, increasing to does not since more information is taken from the existing time 
series. 
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Since the experiment is operating in closed loop, the number of available system identification methods 
is somewhat limited. In this paper we make use of the direct method for closed-loop identification due to 
its simplicity and so identify the unknown parameters of Eq. (15) with least squares, thus minimising the 
sum-of-squares of the model error e(T). However, other methods such as joint input-output identification 
can be used if required [22, 25]. In order to provide sufficiently informative results for system identification 
purposes, the random perturbation r]{t) should have a sufficiently broadband spectrum. However, to minimise 
extraneous noise (that is, random perturbations which are not captured by the discretisation) the bandwidth 
of the random perturbation should be less than the Nyquist frequency corresponding to the discretisation in 
Eq. (15). For the purposes of this paper, ri{t) is generated by passing Gaussian white noise through a 6th 
order Butterworth filter with a cut-off frequency of 10 Hz. 

When there is a significant amount of measurement noise or unmeasurable random disturbances, a 
non-trivial noise model is required in Eq. (15), giving rise to a moving-aver age (MA) term. In this case the 
the unknown coefficients of the resulting ARMAX model must be estimated using a method such as the 
prediction error method (PEM) since the straightforward use of linear least-squares will result in bias [22, 
Chapter 10[. However, linear least-squares provides a quick and effective way of starting the iterative PEM 
optimisation. 

Once a linearised model has been identified the Floquet multipliers of the periodic orbit can be determined 
from the matrix Specifically, we seek to determine the monodromy matrix M such that y{T) = 

Mq~^y(T) with k(T) = 0 and e(r) = 0, that is, we seek a linear mapping which takes one period of data 
points and returns the following period of data points subject to no disturbance to the system input. For the 
method described here M corresponds to the first n rows and n columns of the matrix given by 

R-i(0)(H(l)-i?(0)), (18) 

and the Floquet multipliers of the periodic orbit are the eigenvalues of M. In addition to the Floquet 
multipliers, the eigenvectors of M correspond to the stable and unstable directions of the periodic orbit at a 
particular point in the oscillation. 


4. Experimental apparatus 

To test the effectiveness of this methodology outlined in this paper we apply control-based continuation 
to a nonlinear tuned mass damper (NTMD) similar to the one described in [26[. The NTMD consists of a 
mass able to move horizontally on a low friction bearing system while restrained by two springs that are 
mounted perpendicularly to the direction of motion, thus providing a geometric nonlinearity. The NTMD is 
then excited at the base. This configuration results in a hardening spring-type characteristic. A photograph 
of the experiment is shown in Fig. 2(a) along with a schematic of the experiment in Fig. 2(b). 

The details of the actuation and measurement equipment are as follows. The NTMD is excited using 
an APS-113 long-stroke electrodynamic shaker in current control mode using a Maxon ADS-50/10-4QDC 
motor controller. Typical base displacements are sinusoidal with a frequency ranging from 2.2-3.2 Hz and 
an amplitude ranging from 0-25 mm. The peak response amplitude is limited to ±80 mm; at resonance, 
this limitation restricts the amplitude of the base motion to approximately ±7 mm. The motion of the base 
and the moving mass are measured using laser displacement sensors (an Omron ZX2-LD100 and an Omron 
ZX1-LD300 respectively). In addition to the displacement measurements, the force provided by the shaker is 
measured using an MCL-type load cell. 

For the real-time control, a linear proportional-derivative (PD) controller is used with manually tuned 
gains. The methodology is relatively insensitive to the control gains used provided they are sufficient to 
stabilise any unstable orbits that are encountered. The controller is implemented on a BeagleBone Black 
fitted with a custom data acquisition board (hardware schematics and associated software are open source 
and freely available [27[). All measurements are made at IkHz with no filtering. 

A random perturbation signal is generated, when necessary, on the real-time control board using the 
Box-Muller transformation to generate Gaussian pseudo-random numbers which are then filtered using a 
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Figure 2: Panel (a) shows photograph of the nonlinear tuned mass damper (NTMD) used to test the effectiveness of the 
methodology outlined in this paper. Panel (b) shows a schematic of the NTMD with the springs mounted perpendicular to 
the direction of motion which results in a geometric nonlinearity. Panel (c) shows a schematic of the overall experimental 
rig; the feedback-control loop and a limited amount of signal processing are implemented in real-time while the numerical 
continuation routines are implemented off-line (that is, there are no time constraints on the computations for the continuation; 
the experiment will simply continue running until new input parameters are available). 


sixth-order Butterworth filter with a cut-off arbitrarily set at 10 Hz (below the Nyquist frequency of the 
discretisation used for Eq. (15)). 

Estimations of the Fourier coefficients of the response and the control action are also calculated in 
real-time on the control board. However, this was for convenience rather than a necessity. 


5. Experimental results 

The basic control-based continuation algorithm described in Sec. 2 was used to do repeated continuations 
in the forcing amplitude (the amplitude of displacement of the shaking table) for fixed values of the forcing 
frequency. The forcing frequency, while fixed for individual continuation runs, was varied between 2.2 Hz and 
3.2 Hz in steps of 0.025 Hz. At each data point, full time series measurements were made. These are shown 
as black dots in Fig. 3 where the forcing frequency and forcing amplitude (in mm) are plotted against the 
response amplitude, which we define as the magnitude of the first component in the Fourier series. 

To aid visualisation, a continuous surface constructed from the individual data points is also plotted in 
Fig. 3. This continuous surface is created using Gaussian progress regression on the collected data points 
where the hyper-parameters for the Gaussian process are calculated by maximising the marginal likelihood of 
the hyper-parameters [28, Sec. 5.4]. 

The use of Gaussian process regression (or any other similar scheme for interpolating the multi-dimensional 
data) also allows for geometric features of the solution surface to be easily extracted. One pertinent feature 
is the fold in the solution surface which, from dynamical systems theory, indicates a change in stability of 
the periodic orbits. As such, the fold curve shown in Fig. 3 (blue curve) was extracted using numerical 
continuation with CoCo [20] directly on the regression surface defined by the Gaussian process. 

Other features of interest are frequency response curves which can also be obtained through numerical 
continuation on the regression surface by fixing the forcing amplitude to a prescribed value. Figure 4 shows 
such frequency response curves, including unstable periodic orbits, for fixed forcing amplitudes of F = 1.9, 2, 
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Figure 3: Measurements taken from 41 continuation runs that vary the forcing amplitude for different fixed values of the forcing 
frequency. The forcing frequency is changed in steps between 2.2 Hz and 3.2 Hz. The shaded surface is calculated using Gaussian 
process regression on the measured data points. The location of the unstable periodic orbits in this figure can be inferred from 
the geometry of the solution surface; the plotted curve (solid line) represents a ID fold curve inside which are the unstable 
periodic orbits. This fold curve is calculated using the Gaussian process regressor. 



Figure 4; A series of frequency response curves extracted from the data shown in Fig. 3 using Gaussian process regression. The 
forcing amplitudes of the base motion are 1.9, 2, 2.5, 3 and 3.5 mm respectively. Fold points (limit points) which determine the 
hysteresis region are marked with black dots. 







Figure 5: Two separate time-series measurements in taken open-loop conditions starting from the same unstable periodic orbit. 
The measurements where synchronised using the periodic forcing as a reference signal. Noise in the experiment randomly perturbs 
the trajectory from the unstable periodic orbit to either the stable high-amplitude orbit (blue) or the stable low-amplitude orbit 
(red). This shows that the stable manifold of the unstable periodic orbit acts as a separatrix between the basins of attraction of 
the high- and low-amplitude orbits. 


(a) 10 




Figure 6: Panels (a) and (b) show 2D and 3D state-space projections of a single stable periodic orbit (red) and the randomly 
perturbed orbit (blue) used to calculate the stability of the periodic orbit. Time-delay coordinates are used as a proxy for the 
derivative of the position x{t) to recreate the state space. The time coordinate in panel (b) is normalised such that a period of 
forcing takes one time unit. 


2.5, 3 and 3.5 mm. For high-amplitude forcing, there is little to distinguish the results to those obtained 
from a Duffing equation with hardening nonlinearity. However, for low-amplitude forcing there seems to be a 
significant influence from frictional nonlinearities in the bearing system suspending the mass. 

In order to verify that the unstable orbits found in the experiment are true unstable periodic orbits and 
not artefacts of the control scheme, we repeatedly drive the system to a particular unstable periodic orbit 
and then turn off the stabilising controller. As can be seen from the time series shown in Fig. 5, starting 
from the unstable periodic orbit both the stable low-amplitude and stable high-amplitude periodic orbits can 
be reached — the stable manifold of the unstable orbit acts as a separatrix between the two stable orbits 
as we would expect from dynamical systems theory. Out of 40 separate time series recorded starting from 
the same unstable periodic orbit, 8 end at the high-amplitude periodic orbit and the remaining end at the 
low-amplitude periodic orbit. 

In order to apply the method outlined in Sec. 3, once an orbit has been obtained using control-based 
continuation it must be perturbed with a random input signal. One such example with a perturbation size of 
0.5 is shown in Fig. 6. (Strictly speaking, the perturbation size is measured in volts as it is an input to the 
shaker; however, the spectral content of the perturbation combined with the non-trivial frequency response 
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Figure 7: The mean and 95% confidence interval of the absolute values of the dominant Floquet multipliers estimated against the 
size of the perturbation applied. Each point represents 10 separate measurements. The points marked (a) (red) are calculated 
for a stable periodic orbit in open-loop conditions; the points marked (b) (blue) are calculated for the same stable periodic 
orbit in closed-loop conditions; and the points marked (c) (green) are calculated for an unstable periodic orbit in closed-loop 
conditions. For the stable periodic orbits (a) and (b) the Floquet multipliers are complex; the estimated errors in the real and 
imaginary parts are approximately equal. 


of the shaker mean that it is not straightforward or useful to state the perturbation size in mm.) To avoid 
estimating the velocity of motion, a time-delay coordinate x{t — t) is used to reconstruct the state-space of 
the experiment. Here the value of r used is T/5 where T is the period of the forcing. 

There is a single parameter in Sec. 3 for which there is no algorithmic way to determine an appropriate 
value, that is the amplitude of the perturbation applied to the periodic orbit. Thus in order to determine an 
appropriate amplitude we select two periodic orbits, one stable and one unstable and calculate the Floquet 
multipliers using the fitted linear time-varying model Eq. (15) for a variety of perturbation sizes and repeat 
the experiments 10 times. 

Throughout this paper we set m = 10 and n = 4 in Eq. (15). 

The results of the Floquet multiplier estimations are shown in Fig. 7. The points marked (a) (in red) are 
the absolute values of the Floquet multipliers estimated for a stable periodic orbit while the experiment is 
running in open-loop; similarly, the points marked (b) (in blue) are the Floquet multipliers estimated for the 
same stable periodic orbit while the experiment is running in closed-loop. Finally, the points marked (c) (in 
green) are the Floquet multipliers estimated for an unstable periodic orbit while the experiment is running 
in closed-loop. Each Floquet multiplier is marked with the 95% confidence range. 

For both the stable and the unstable periodic orbits it can be seen that the 95% confidence range 
narrows significantly for larger perturbation sizes; the larger perturbations allow the Floquet multipliers to 
be estimated more consistently. However, the consistency of the results does not imply accuracy — as with 
estimating derivatives from finite differences, taking a large step introduces errors caused by the nonlinearities 
in the system. Though, since the mean magnitude of the Floquet multipliers does not change significantly 
beyond a perturbation size of 0.2, we can have reasonable confidence in the results and so a perturbation size 
of 0.5 is used throughout the remainder of this paper. 

Furthermore, Fig. 7 shows that running the system in closed-loop rather than open-loop does not have 
a significant affect on the estimation of the Floquet multipliers — the error bars of the points (a) and (b) 
overlap considerably. 

Unfortunately no accurate independent estimations of the Floquet multipliers are available to check the 
results of the system identification. Convergence tests on the stable periodic orbit were performed, however 
the transient dynamics of the electrodynamic shaker as the experiment equilibrated rendered the results 
meaningless. Furthermore, escape tests from the unstable periodic orbit such as those seen in Fig. 5 produced 
estimates ranging from 1.2 to 1.5 depending on the measurement cut-offs used. As such, we consider the 
best way to judge the accuracy of the estimations is via comparison with geometric phenomena such as fold 
points where a Floquet multiplier should pass through -|-1 in the complex plane. 

Subsequently, a continuation in the forcing amplitude was performed with Floquet multipliers estimated 


10 










2 4 6 8 

Forcing amplitude (mm) 



Figure 8: Panel (a) shows a single continuation with the response amplitude plotted against the forcing amplitude for a fixed 
fixed forcing frequency of 3.2 Hz. Stable periodic orbits are shown by solid blue dots and unstable orbits are shown by red circles. 
The maximum displacement of the experiment is limited to ±80 mm. Panel (b) shows the absolute value of the estimated 
Floquet multipliers plotted against the response amplitude for three separate continuation runs (each shown with a different 
symbol). It can be seen that most of the stable periodic orbits have complex conjugate Floquet multipliers which then become 
real close to the saddle-node bifurcations; between the two saddle-node bifurcations, one of these Floquet multipliers lies outside 
the unit circle while the other lies inside the unit circle. 


for each obtained solution; the results are shown in Fig. 8(a) with stable periodic orbits marked as solid blue 
dots and unstable periodic orbits marked as red circles. The results agree very well with what is expected 
from dynamical systems theory; the orbits between the two fold points are all unstable. 

To ensure repeatability, the same continuation run was carried out three times and the absolute values of 
the corresponding Floquet multipliers are plotted in Fig. 8(b) with different symbols for the different runs. 
As can be seen from Fig. 8(b) the results are very consistent with all the runs showing good quantitative as 
well as qualitative agreement. As expected, a complex conjugate pair of Floquet multipliers becomes real 
close to the fold point after which a single real multiplier passes through the unit circle. This process reverses 
close to the second fold point. 

As well as providing information about the Floquet multipliers of the periodic orbits, the monodromy 
matrix obtained from Eq. (15) also provides the stable and unstable eigendirections of the periodic orbit. 
Figure 9 shows the time series from Fig. 5 plotted in state-space, with the unstable periodic orbit (marked in 
red) lying between the two stable periodic orbits (marked in green). From Fig. 9(b) it can be seen that the 
escape from the unstable periodic orbit occurs along a two-dimensional unstable manifold. 

Figure 10 shows a Poincare section at a fixed time in the forcing cycle that was created with the data 
from Fig. 9 combined with data from a further 38 independent time series measurements; the blue dots 
denote intersections of the trajectories with the Poincare section. The intersection of the two-dimensional 
unstable manifold with the Poincare section results in a well defined one-dimensional curve showing how 
trajectories leave the unstable periodic orbit (marked with a red dot) and approach the two stable periodic 
orbits (marked with green crosses). 

Superimposed on Fig. 10 is a set of red arrows which mark the unstable and stable eigendirections 
calculated from the monodromy matrix of the periodic orbit using only local perturbations (as previously 
described). The eigendirections represent the linearisation of the manifolds at the unstable periodic orbit 
and, as such, show remarkably good agreement with the measured unstable manifold which was calculated in 
open-loop conditions. This opens up the possibility for using eigendirection information in other calculations 
on the physical system, for example the estimation of basins of attraction. 


6. Conclusions 

This paper has proposed a new method for estimating Floquet multipliers and their associated eigendirec¬ 
tions for periodic orbits that encountered when using control-based continuation on a physical experiment. A 
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Figure 9: State-space projections in 2D and 3D of the data shown in Fig. 5 using time-delay coordinates. The red curve 
represents the unstable periodic orbit and the two green curves represent the stable high- and low-amplitude periodic orbits. 
The blue curves show the transient dynamics between the unstable orbits and the stable orbits. The two different stable orbits 
are reached with approximately equal probability. The time coordinate in panel (b) is normalised such that a period of forcing 
takes one time unit. 



Figure 10: A Poincare section showing measurements from 40 separate time series of the escape from the unstable periodic orbit 
shown in Fig. 5 and Fig. 9. Each dot (blue) corresponds to an intersection with the Poincare section defined by t mod T = 0, 
where T is the forcing period. The dots trace out the ID unstable manifold of the unstable orbit. The positions of the stable 
high- and low-amplitude periodic orbits are marked with green crosses. The stable and unstable eigenspaces estimated using the 
methodology outlined in this paper are marked by red arrows. The unstable eigenspace shows good agreement with the dots 
from the time series measurements. 
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linear time-varying model in the form of an auto-regressive model with exogenous inputs (ARX) is fitted to 
each periodic orbit using small perturbations to the orbit to explore the nearby state-space. 

The method was demonstrated on a nonlinear mass-spring-damper-type experiment that has a hardening 
spring characteristic. The nonlinearity is provided by placing springs perpendicular to the direction of motion, 
thus creating a geometric nonlinearity. The Floquet multiplier estimations were shown to agree with what 
is expected from dynamical systems theory and the associated eigendirections match well with open-loop 
measurements taken. 
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All the experimental data used in this paper has been deposited into the University of Bristol Research 
Data Repository and is publicly available for download [29]. 
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